* Figure4.do  09/05/22
* blacks and whites change ranks of opioid death rate

do Figure4prep  // Note that, as CDC does, 2021 pop is copied from 2020
keep if ucd=="Drug" & mcd=="Opioid" & year>=2012 & (age=="0-44" | age=="45-64" | age=="65+")
gen opioiddeathrate = deaths*100000/pop
gen byte black = race==`"{"Black or African American"}"'
gen byte female = gender==`"{"Female"}"'
encode age, gen(agenum)
encode state, gen(censusdiv)
gen gapraw = .
gen gapadj = .
gen gapold = .
gen gapyoung = .
gen gapmid = .
gen gapeld = .
gen seadj = .
set more off
foreach yr of numlist 2012/2021 {
  reg opioiddeathrate black if year==`yr' [fw=pop]
  replace gapraw = _b[black] if year==`yr'
  reg opioiddeathrate black female i.agenum i.censusdiv if year==`yr' [fw=pop]
  replace gapadj = _b[black] if year==`yr'
  replace seadj = _se[black] if year==`yr'
  reg opioiddeathrate black female i.agenum i.censusdiv if year==`yr' & age~="0-44" [fw=pop]
  replace gapold = _b[black] if year==`yr'
    reg opioiddeathrate black female i.censusdiv if year==`yr' & age=="0-44" [fw=pop]
  replace gapyoung = _b[black] if year==`yr'
  reg opioiddeathrate black female i.censusdiv if year==`yr' & age=="45-64" [fw=pop]
  replace gapmid = _b[black] if year==`yr'
  reg opioiddeathrate black female i.censusdiv if year==`yr' & age=="65+" [fw=pop]
  replace gapeld = _b[black] if year==`yr'
}
gen zeroes = 0
#delimit ;
twoway (line zeroes year if black==0 & female==0 & agenum==1 & censusdiv==1, lcolor(gray) lpattern(dash) xlabel(2012[1]2021))
  (connected gapadj gapold gapyoung year if black==0 & female==0 & agenum==1 & censusdiv==1, lcolor(black green blue) mcolor(black green blue) msymbol(circle square triangle) msize(medium medium medium)
  title("Figure 4.  Black-white gaps in opioid fatality rates") subtitle("adjusted for gender, age group, and Census division")
  legend(order(3 "Ages 45+" 2 "All ages" 4 "Ages 0-44")  region(lcolor(white)))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Deaths per 100K population")
  note("Notes: For each year an opioid death indicator is regressed on indicators for race, gender, age group, and Census division, with"
  "100K times the black coefficient shown in the figure.  A regression observation is a black or white U.S. resident alive January 1."
  "Gap confidence intervals (not shown) are less then 2 per 100K population.  Population and mortality source: CDC Wonder.", size(vsmall)) );
#delimit cr


* fw=pop confidence intervals are miniscule.  CIs are larger when estimated with individual-level data (each year having one record for each person alive at the beginning of the year)
* see Figure4stderr.do, which calculates the CIs cited in the twoway note()

tabstat pop deaths if black==1, by(year) stat(sum) format(%12.0fc)

